#include "numb.h"

using namespace std;
/*
 * @brief Barrett模约减流程，根据相关预计算量计算  a mod q
 * 
 * @param a 要约减的数a，长度为128bit
 * @param q 模数q,长度64bit
 * @param muHigh Barrett约减中的预计算量mu的高64位，mu为2^128−1除以n的商
 * @param muLow Barrett约减中的预计算量mu的低64位
 * @return uint64_t a mod q
 */
uint64_t SingleBarrettNew(__uint128_t& a, uint64_t q, uint64_t muHigh, uint64_t muLow)
{
    return Barrett128(a, q, muHigh, muLow);
}

/**
 * 功能描述：Barrett计算流程的汇编优化，供上层Barrett函数调用
 * @__uint128_t& a:要约减的数uint128bit的数a
 * @uint64_t q：模数q
 * @uint64_t mu_high：Barrett约减中的预计算量mu的高64位，mu为2^128−1除以n的商
 * @uint64_t mu_low：Barrett约减中的预计算量mu的低64位
 * 
*/
uint64_t Barrett128(__uint128_t &a, uint64_t q, uint64_t mu_high, uint64_t mu_low)
{
    uint64_t a_high = a >> 64;
    uint64_t a_low = static_cast<uint64_t>(a);
    uint64_t qq = 2 * q;
    asm volatile (
        // Step 1: Calculate a_low * mu_high, store high 64 bits in x1
            "umulh x1, %[a_low], %[mu_high]\n\t"

            // Step 2: Calculate a_high * mu_low and add to x1
            "umulh x0, %[a_high], %[mu_low]\n\t"
            "add x1, x1, x0\n\t"

            // Step 3: Calculate a_high * mu_high and add to x1
            "mul x0, %[a_high], %[mu_high]\n\t"
            "add x1, x1, x0\n\t"

            // Step 4: Calculate x1 * q
            "mul x1, x1, %[q]\n\t"

            // Step 5: Calculate a_low - x1
            "sub x1, %[a_low], x1\n\t"

            // Step 6: if (x1 >= qq) x1 -= qq
            "subs x0, x1, %[qq]\n\t"  // Calculate x1 - qq and set flags
            "csel x1, x0, x1, hs\n\t"  // If x1 >= qq (hs), select x1 - qq, else select x1

            // Step 7: if (x1 >= q) x1 -= q
            "subs x0, x1, %[q]\n\t"  // Calculate x1 - q and set flags
            "csel %[a_low], x0, x1, hs\n\t"  // If x1 >= q (hs), select x1 - q, else select x1

        // Store the result back to a_low
        // "mov %[a_low], x1\n\t"

            : [a_low]"+r"(a_low)  // Output: a_low (modified)
    : [a_high]"r"(a_high), [mu_low]"r"(mu_low), [mu_high]"r"(mu_high), [q]"r"(q), [qq]"r"(qq)
    : "x0", "x1", "memory", "cc"
    );


    return a_low;
}

/**
 * 功能描述：对 64-bit 的 a 进行 Barrett 约减
 * @uint64_t a: 待约减的数
 * @uint64_t q：模数 q
 * @uint64_t mu：Barrett 约减的预计算量，mu = floor(2^64 / q)
 */
uint64_t Barrett64(uint64_t a, uint64_t q, uint64_t mu)
{
    asm volatile (
        // Step 1: Calculate a * mu, store high 64 bits in x1
            "umulh x1, %[a], %[mu]\n\t"

            // Step 2: Calculate x1 * q
            "mul x1, x1, %[q]\n\t"

            // Step 3: Calculate a - x1
            "sub x1, %[a], x1\n\t"

            // Step 4: if (x1 >= q) x1 -= q
            "subs x0, x1, %[q]\n\t"  // Calculate x1 - q and set flags
            "csel %[a], x0, x1, hs\n\t"  // If x1 >= q (hs), select x1 - q, else select x1

            : [a]"+r"(a)  // Output: a (modified)
    : [q]"r"(q), [mu]"r"(mu)  // Inputs: q, mu
    : "x0", "x1", "cc"  // Clobbers: x0, x1, condition flags
    );
    return a;
}

/**
 * 功能描述：对 64-bit 的 a 进行 Barrett 约减
 * @uint64x2_t a：存在NEON寄存器中的待约简的两个64bit数a
 * @uint64x2_t q：存在NEON寄存器中的待约简的两个模数q
 * @uint64x2_t mu：Barrett 约减的预计算量，mu
 * @uint64x2_t *res: 两个数的约简结果
 */
void Barrett64NEON(uint64x2_t *res, uint64x2_t a, uint64x2_t q, uint64x2_t mu)
{
    uint64_t res0 = 0, res1 = 0;

    asm volatile (
            "umulh x1, %[a0], %[mu0]\n\t"
            "mul x1, x1, %[q0]\n\t"
            "sub x1, %[a0], x1\n\t"
            "subs x0, x1, %[q0]\n\t"
            "csel %[res0], x0, x1, hs\n\t"

            : [res0]"=r"(res0)
    : [a0]"r"(vgetq_lane_u64(a, 0)), [q0]"r"(vgetq_lane_u64(q, 0)), [mu0]"r"(vgetq_lane_u64(mu, 0))
    : "x0", "x1", "cc"
    );

    asm volatile (
            "umulh x1, %[a1], %[mu1]\n\t"
            "mul x1, x1, %[q1]\n\t"
            "sub x1, %[a1], x1\n\t"
            "subs x0, x1, %[q1]\n\t"
            "csel %[res1], x0, x1, hs\n\t"

            : [res1]"=r"(res1)
    : [a1]"r"(vgetq_lane_u64(a, 1)), [q1]"r"(vgetq_lane_u64(q, 1)), [mu1]"r"(vgetq_lane_u64(mu, 1))
    : "x0", "x1", "cc"
    );

    *res = vsetq_lane_u64(res0, *res, 0);
    *res = vsetq_lane_u64(res1, *res, 1);

}
/**
 * @uint64_t a：需要求相反数的值
 * @uint64_t& r：求相反数的结果
 */
void Negate(uint64_t& r, uint64_t a)
{
    r = -a;
}

/**
 * @uint64_t& r:返回结果(a+b)%p
 * @uint64_t a:模加的操作数a
 * @uint64_t b：模加的操作数b
 * @uint64_t m：模数p
*/
void AddMod(uint64_t& r, uint64_t a, uint64_t b, uint64_t p)
{
    ComSubQ(a, p);
    ComSubQ(b, p);
    r = a + b;
    if (r >= p){
        r -=p;
    }
}
/**
 * @uint64_t& r:返回结果(a - b)%p
 * @uint64_t a:模减的操作数a
 * @uint64_t b：模减的操作数b
 * @uint64_t p：模数p
*/
void SubMod(uint64_t& r, uint64_t a, uint64_t b, uint64_t p)
{
    if (p == 0) {
        throw std::runtime_error("Divisor mod by zero is not allowed.");
    }
    ComSubQ(a, p);
    ComSubQ(b, p);
    r = b % p;
    r = (a + p - r) % p;
}
/*
 * @uint64_t& r:返回结果(a * b)%p
 * @uint64_t a:模乘的操作数a
 * @uint64_t b：模乘的操作数b
 * @uint64_t p：模数p
 */
void MulMod(uint64_t& r, uint64_t a, uint64_t b, uint64_t p)
{
    if (p == 0) {
        throw std::runtime_error("Divisor mod by zero is not allowed.");
    }
    if (a > p) {
        a = a % p;
    }
    if (b > p) {
        b = b % p;
    }
    uint128_t mul = static_cast<uint128_t>(a) * (b);
    mul %= static_cast<uint128_t>(p);
    r = static_cast<uint64_t>(mul);
}

/*
 * 功能描述：计算 a*b mod p，并且通过barrett约减来计算最终结果
 * @uint64_t& r：用来传递计算结果
 * @uint64_t a: 源操作数a
 * @uint64_t b: 源操作数b
 * @uint64_t p: 模数p
 * @uint64_t twop：模数p的两倍
 * @uint64_t prHigh: Barrett预计算量Pr的高64bit       #在Context中定义对应的Pr矩阵
 * @uint64_t prLow:  Barrett预计算量Pr的低64bit
 */
void MulModBarrett(uint64_t& r, uint64_t a, uint64_t b, uint64_t p, uint64_t twop, uint64_t prHigh, uint64_t prLow)
{
    if (p == 0) {
        throw std::runtime_error("Divisor mod by zero is not allowed.");
    }
    if (a > p) {
        a = a % p;
    }
    if (b > p) {
        b = b % p;
    }
    asm volatile (
        //计算 a*b  低位存x3 高位存x3寄存器
            "mul   x3 ,%[a] , %[b] \n\t "
            "umulh x4 ,%[a] , %[b] \n\t"
            //计算Mul_l 与预计算量pr的乘积，高64位存在res中
            "umulh x1, x3, %[mu_high]\n\t"

            // 计算 a.high * mu.low + res : 先计算高64位乘法，再累加res,要存放到res中
            "umulh x0, x4, %[mu_low]\n\t"
            "add x1, x1, x0\n\t"

            // 计算 a.high * mu.high + res,先计算低64位乘法， 再累加res,要存放到res中
            "mul x0, x4, %[mu_high]\n\t"
            "add x1, x1, x0\n\t"

            // 计算 res * q 存回res
            "mul x1, x1, %[q]\n\t"

            // 计算 a.low - res,存回res
            "sub x1, x3, x1\n\t"

            // 计算res 和2*q的比较，如果直接减发生溢出就set flag
            "subs x0, x1, %[qq]\n\t"
            "csel x1, x0, x1, hs\n\t"

            // 计算res 和q的比较，如果直接减发生溢出就set flag
            "subs x0, x1, %[q]\n\t"
            "csel %[r], x0, x1, hs\n\t"
        // Store the result back to a_low
        //   "mov %[r], x1\n\t"
            : [a]"+r"(a) , [b]"+r"(b) ,[r]"+r"(r)
    :  [mu_low]"r"(prLow), [mu_high]"r"(prHigh), [q]"r"(p), [qq]"r"(twop)
    : "x0", "x1","x3","x4", "memory", "cc"
    );
}
/*
 * 功能描述：计算 a 的barrett约减来计算最终结果
 * @uint64_t& r：用来传递计算结果
 * @uint64_t a: 源约简操作数a
 * @uint64_t p: 模数m
 * @uint64_t twop：模数p的两倍
 * @uint64_t prHigh: Barrett预计算量Pr的高64bit       #在Context中定义对应的Pr矩阵
 * @uint64_t prLow:  Barrett预计算量Pr的低64bit
 */
void ModBarrett(uint64_t &r, uint64_t a, uint64_t p, uint64_t prHigh, uint64_t prLow)
{
    __uint128_t a_ = static_cast<__uint128_t>(a);
    r = Barrett128(a_ ,p, prHigh, prLow);
}
void ModBarrett(uint64_t &r, uint128_t a, uint64_t p, uint64_t prHigh, uint64_t prLow)
{
    r = Barrett128(a, p, prHigh, prLow);
}
/*
 * @uint64_t x：模逆的操作数x
 * @uint64_t p：模数p
 * @return  :返回结果x^-1%p
 */
uint64_t InvMod(uint64_t x, uint64_t p)
{
    int64_t m0 = p;
    int64_t x0 = 0, x1 = 1;
    int64_t q, t;

    if (p == 1)
        return 0; // 任何数 mod 1 都没有逆元

    while (x > 1) {
        // q 是商
        q = x / p;
        t = p;
        // m 是余数，a 是被除数
        p = x % p;
        x = t;
        // 更新 x0 和 x1
        t = x0;
        x0 = x1 - q * x0;
        x1 = t;
    }

    // 确保 x1 是正数
    if (x1 < 0)
        x1 += m0;

    if (x1 == -1) {
        //        throw invalid_argument("Inverse doesn't exist!!!");
    }

    return x1;
}
/*
 * @brief 计算x^y mod p
 * @param x
 * @param y
 * @param p
 * @return uint64_t
 */
uint64_t PowMod(uint64_t x, uint64_t y, uint64_t p)
{
    uint64_t res = 1;
    while (y > 0) {
        if (y & 1) {
            MulMod(res, res, x, p);
        }
        y = y >> 1;
        MulMod(x, x, x, p);
    }
    return res;
}

/*
 * @brief ：计算 x 的乘法逆元 (mod 2^64)
 * 
 * @param x ：输入值
 * @return uint64_t ：返回 x 的乘法逆元
 */
uint64_t Inv(uint64_t x)
{
    return Pow(x, static_cast<uint64_t>(-1));
}

/*
 * @brief 计算 x 的 y 次幂 (x^y)
 * @param x 底数
 * @param y 指数
 * @return uint64_t 返回 x^y 的结果
 */
uint64_t Pow(uint64_t x, uint64_t y)
{
    uint64_t res = 1;
    while (y > 0) {
        if (y & 1) {
            res *= x;
        }
        y = y >> 1;
        x *= x;
    }
    return res;
}

/*
 * @brief 计算 32 位无符号整数的位反转（Bit Reverse）
 * @param x 输入的 32 位无符号整数
 * @return uint32_t 返回位反转后的值
 */
uint32_t BitReverse(uint32_t x)
{
    x = (((x & 0xaaaaaaaa) >> 1) | ((x & 0x55555555) << 1));
    x = (((x & 0xcccccccc) >> 2) | ((x & 0x33333333) << 2));
    x = (((x & 0xf0f0f0f0) >> 4) | ((x & 0x0f0f0f0f) << 4));
    x = (((x & 0xff00ff00) >> 8) | ((x & 0x00ff00ff) << 8));
    return ((x >> 16) | (x << 16));
}

/*
 * @brief 计算两个无符号整数的最大公约数 (GCD)
 * @param a 第一个无符号整数
 * @param b 第二个无符号整数
 * @return uint64_t 返回 a 和 b 的最大公约数
 */
uint64_t Gcd(uint64_t a, uint64_t b)
{
    if (a == 0) {
        return b;
    }
    return Gcd(b % a, a);
}

/*
 * @brief 计算两个有符号整数的最大公约数 (GCD)
 * @param a 第一个有符号整数
 * @param b 第二个有符号整数
 * @return long 返回 a 和 b 的最大公约数
 */
long Gcd(long a, long b)
{
    if (a == 0) {
        return b;
    }
    return Gcd(b % a, a);
}

/*
 * @brief 查找一个数的所有素因子
 * @param s 用于存储素因子的集合
 * @param number 输入的整数
 */
void FindPrimeFactors(set<uint64_t>& s, uint64_t number)
{
    while (number % 2 == 0) {
        s.insert(2);
        number /= 2;
    }
    for (uint64_t i = 3; i < sqrt(number); i++) {
        while (number % i == 0) {
            s.insert(i);
            number /= i;
        }
    }
    if (number > 2) {
        s.insert(number);
    }
}

/*
 * @brief 查找模数的一个原根
 * @param p 模数
 * @return uint64_t 返回找到的原根，如果不存在则返回 -1
 */
uint64_t FindPrimitiveRoot(uint64_t p)
{
    set<uint64_t> s;
    uint64_t phi = p - 1;
    FindPrimeFactors(s, phi);
    for (uint64_t r = 2; r <= phi; r++) {
        bool flag = false;
        for (auto it = s.begin(); it != s.end(); it++) {
            if (PowMod(r, phi / (*it), p) == 1) {
                flag = true;
                break;
            }
        }
        if (flag == false) {
            return r;
        }
    }
    return -1;
}

/*
 * @brief 查找模数的第 M 个单位根
 * @param M 单位根的阶数
 * @param p 模数
 * @return uint64_t 返回找到的单位根，如果不存在则返回 -1
 */
uint64_t FindMthRootOfUnity(uint64_t M, uint64_t p)
{
    if (M == 0) {
        throw std::runtime_error("Divisor mod by zero is not allowed.");
    }
    uint64_t res;
    res = FindPrimitiveRoot(p);
    if ((p - 1) % M == 0) {
        uint64_t factor = (p - 1) / M;
        res = PowMod(res, factor, p);
        return res;
    }
    else {
        return -1;
    }
}

/*
 * @brief 使用 Shoup 模乘的预计算参数
 * @param x 输入值
 * @param p 模数
 * @return uint64_t 返回计算结果
 */
uint64_t xShoup(uint64_t x, uint64_t p)
{
    if (p == 0) {
        throw std::runtime_error("Divisor mod by zero is not allowed.");
    }
    uint128_t temp = x;
    temp <<= 64;
    return uint64_t(temp / p);
}

/*
 * @brief 使用 Shoup 模乘计算乘法
 * @param a 第一个乘数
 * @param b 第二个乘数
 * @param b_shoup 第二个乘数的预计算值
 * @param q 模数
 * @return uint64_t 返回模乘结果
 */
uint64_t ShoupMul(uint64_t a, uint64_t b, uint64_t b_shoup, uint64_t q)
{
    uint128_t v_w_mul = static_cast<uint128_t>(a) * b_shoup;
    uint64_t hi = v_w_mul >> 64;
    uint64_t res = a * b - hi * q;
    ComSubQ(res, q);

    // uint64_t res0 = 0;
    // MulMod(res0, a, b, q);
    // if(res0 != res) {
    //     printf("!=\n");
    // }

    return res;
}

/*
 * @brief 使用 Shoup 模乘计算乘法
 * @param res 用于存储结果的引用
 * @param a 第一个乘数
 * @param b 第二个乘数
 * @param bShoup 第二个乘数的预计算值
 * @param q 模数
 */
void MulModShoup(uint64_t &res, uint64_t a, uint64_t b, uint64_t bShoup, uint64_t q)
{
    uint128_t v_w_mul = static_cast<uint128_t>(a) * bShoup;
    uint64_t hi = v_w_mul >> 64;
    res = a * b - hi * q;
    if (res >= q)
        res -= q;
}

/*
 * @brief 使用 Barrett Reduction 和 NEON 指令计算模乘
 * @param r 用于存储结果的 NEON 向量指针
 * @param a 第一个乘数
 * @param b 第二个乘数
 * @param p 模数
 * @param twop 2 * p，用于快速比较
 * @param prHigh 预计算量的高位部分
 * @param prLow 预计算量的低位部分
 */
void MulModBarrettNeon(uint64x2_t* r, uint64x2_t a, uint64x2_t b, uint64x2_t p, uint64x2_t twop, uint64x2_t prHigh, uint64x2_t prLow) {
    uint64_t res0 = 0, res1 = 0;

    // 处理第一个元素
    asm volatile(
        // 计算 a[0] * b[0] 低位存 x3, 高位存 x4
            "mul x3, %[a0], %[b0]\n\t"
            "umulh x4, %[a0], %[b0]\n\t"

            // 计算 Mul_l 与预计算量 prHigh 的乘积，高64位存 x1
            "umulh x1, x3, %[prHigh0]\n\t"

            // 计算 a[0].high * prLow + res (高位乘法再累加到 x1)
            "umulh x0, x4, %[prLow0]\n\t"
            "add x1, x1, x0\n\t"

            // 计算 a[0].high * prHigh + res (低位乘法再累加到 x1)
            "mul x0, x4, %[prHigh0]\n\t"
            "add x1, x1, x0\n\t"

            // 计算 x1 * p 存回 x1
            "mul x1, x1, %[p0]\n\t"

            // 计算 a[0].low - x1, 存回 x1
            "sub x1, x3, x1\n\t"

            // 计算 x1 和 2 * p 的比较，如果直接减发生溢出就 set flag
            "subs x0, x1, %[twop0]\n\t"
            "csel x1, x0, x1, hs\n\t"

            // 计算 x1 和 p 的比较，如果直接减发生溢出就 set flag
            "subs x0, x1, %[p0]\n\t"
            "csel %[res0], x0, x1, hs\n\t"

            : [res0] "=r" (res0)
    : [a0] "r" (vgetq_lane_u64(a, 0)), [b0] "r" (vgetq_lane_u64(b, 0)),
    [p0] "r" (vgetq_lane_u64(p, 0)), [twop0] "r" (vgetq_lane_u64(twop, 0)),
    [prHigh0] "r" (vgetq_lane_u64(prHigh, 0)), [prLow0] "r" (vgetq_lane_u64(prLow, 0))
    : "x0", "x1", "x3", "x4", "cc"
    );

    // 处理第二个元素
    asm volatile(
        // 计算 a[1] * b[1] 低位存 x3, 高位存 x4
            "mul x3, %[a1], %[b1]\n\t"
            "umulh x4, %[a1], %[b1]\n\t"

            // 计算 Mul_l 与预计算量 prHigh 的乘积，高64位存 x1
            "umulh x1, x3, %[prHigh1]\n\t"

            // 计算 a[1].high * prLow + res (高位乘法再累加到 x1)
            "umulh x0, x4, %[prLow1]\n\t"
            "add x1, x1, x0\n\t"

            // 计算 a[1].high * prHigh + res (低位乘法再累加到 x1)
            "mul x0, x4, %[prHigh1]\n\t"
            "add x1, x1, x0\n\t"

            // 计算 x1 * p 存回 x1
            "mul x1, x1, %[p1]\n\t"

            // 计算 a[1].low - x1, 存回 x1
            "sub x1, x3, x1\n\t"

            // 计算 x1 和 2 * p 的比较，如果直接减发生溢出就 set flag
            "subs x0, x1, %[twop1]\n\t"
            "csel x1, x0, x1, hs\n\t"

            // 计算 x1 和 p 的比较，如果直接减发生溢出就 set flag
            "subs x0, x1, %[p1]\n\t"
            "csel %[res1], x0, x1, hs\n\t"

            : [res1] "=r" (res1)
    : [a1] "r" (vgetq_lane_u64(a, 1)), [b1] "r" (vgetq_lane_u64(b, 1)),
    [p1] "r" (vgetq_lane_u64(p, 1)), [twop1] "r" (vgetq_lane_u64(twop, 1)),
    [prHigh1] "r" (vgetq_lane_u64(prHigh, 1)), [prLow1] "r" (vgetq_lane_u64(prLow, 1))
    : "x0", "x1", "x3", "x4", "cc"
    );

    // 将结果存储回 r 向量
    *r = vsetq_lane_u64(res0, *r, 0);
    *r = vsetq_lane_u64(res1, *r, 1);
}

/*
 * @brief 使用 NEON 指令计算模加
 * @param r 用于存储结果的 NEON 向量指针
 * @param a 第一个加数
 * @param b 第二个加数
 * @param p 模数
 */
void AddModNEON(uint64x2_t* r, uint64x2_t a, uint64x2_t b, uint64x2_t p) {
    uint64x2_t sum = vaddq_u64(a, b);
    uint64x2_t cmp = vcgtq_u64(sum, p);
    *r = vsubq_u64(sum, vandq_u64(cmp, p));
}

/*
 * @brief 使用 Shoup 技术和 NEON 指令计算模乘
 * @param result 用于存储结果的 NEON 向量指针
 * @param a 第一个乘数
 * @param b 第二个乘数
 * @param bShoup 第二个乘数的预计算值
 * @param q 模数
 */
void ShoupMulNEON(uint64x2_t* result, uint64x2_t a, uint64x2_t b, uint64x2_t bShoup, uint64x2_t q) {
    uint64_t res0 = 0, res1 = 0;

    // 处理第一个元素
    asm volatile(
            "mul    x3, %[a0], %[bShoup0] \n\t" // x3 = a[0] * bShoup[0] (low 64 bits)
            "umulh  x4, %[a0], %[bShoup0] \n\t" // x4 = a[0] * bShoup[0] (high 64 bits)
            "mul    x1, x4, %[q0]        \n\t" // x1 = (a[0] * bShoup[0] high) * q[0]
            "mul    x0, %[a0], %[b0]     \n\t" // x0 = a[0] * b[0] (low 64 bits)
            "sub    x1, x0, x1           \n\t" // x1 = (a[0] * b[0]) - ((a[0] * bShoup[0] high) * q[0])
            "subs   x0, x1, %[q0]        \n\t" // x0 = x1 - q[0]
            "csel   %[res0], x1, x0, lo  \n\t" // res0 = (x1 < q[0]) ? x1 : x0
            : [res0] "=r" (res0)
    : [a0] "r" (vgetq_lane_u64(a, 0)), [b0] "r" (vgetq_lane_u64(b, 0)), [bShoup0] "r" (vgetq_lane_u64(bShoup, 0)), [q0] "r" (vgetq_lane_u64(q, 0))
    : "x0", "x1", "x3", "x4", "cc"
    );

    // 处理第二个元素
    asm volatile(
            "mul    x3, %[a1], %[bShoup1] \n\t" // x3 = a[1] * bShoup[1] (low 64 bits)
            "umulh  x4, %[a1], %[bShoup1] \n\t" // x4 = a[1] * bShoup[1] (high 64 bits)
            "mul    x1, x4, %[q1]        \n\t" // x1 = (a[1] * bShoup[1] high) * q[1]
            "mul    x0, %[a1], %[b1]     \n\t" // x0 = a[1] * b[1] (low 64 bits)
            "sub    x1, x0, x1           \n\t" // x1 = (a[1] * b[1]) - ((a[1] * bShoup[1] high) * q[1])
            "subs   x0, x1, %[q1]        \n\t" // x0 = x1 - q[1]
            "csel   %[res1], x1, x0, lo  \n\t" // res1 = (x1 < q[1]) ? x1 : x0
            : [res1] "=r" (res1)
    : [a1] "r" (vgetq_lane_u64(a, 1)), [b1] "r" (vgetq_lane_u64(b, 1)), [bShoup1] "r" (vgetq_lane_u64(bShoup, 1)), [q1] "r" (vgetq_lane_u64(q, 1))
    : "x0", "x1", "x3", "x4", "cc"
    );

    // 将结果存储回 r 向量
    *result = vsetq_lane_u64(res0, *result, 0);
    *result = vsetq_lane_u64(res1, *result, 1);

}

/*
 * @brief 计算 a = a - q，如果 a >= q
 * @param a 被减数（引用）
 * @param q 减数
 */
void ComSubQ(uint64_t& a, uint64_t q) {
    if(a >= q) {
        a -= q;
    }
}

/*
 * @brief 使用 NEON 指令计算 a = a - q，如果 a >= q
 * @param a 被减数（NEON 向量引用）
 * @param q 减数（NEON 向量）
 */
void ComSubQNeon(uint64x2_t &a, uint64x2_t q) {
    uint64x2_t mask = vcgeq_u64(a, q); // mask = (a >= q) ? 0xFFFFFFFFFFFFFFFF : 0x0
    uint64x2_t sub_value = vandq_u64(mask, q); // mask AND q (preserve q where condition met)
    a = vsubq_u64(a, sub_value); // a = a - sub_value (only subtract where condition met)
}